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By means of molecular-dynamics simulations, temperature driven diffusionless structural phase 
transitions in equi- and nearly equiatomic ordered nickel-titanium alloys were investigated. For this 
purpose, a model potential from the literature was adopted [W. S. Lai and B. X. Liu, J. Phys. 
Condens. Matter 12, L53 (2000)], which is based on the tight-binding model in second moment 
approximation. The model predicts a stable B19' phase at low temperatures and a nearly cubic B2 
phase at high temperatures. After an analysis of crystallography and energetics of the emerging 
structures, the experimentally known strong dependence of transition temperatures on composition 
is confirmed and related to lattice instability. Free energy calculations finally give insight into the 
driving forces of the phase transitions, and reveal free energy barriers inhibiting them below the 
O ■ transition temperatures. 
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PACS numbers: 02.70.Ns, 05.70.Fh, 64.60.Dc, 81.30.Kf 



I. INTRODUCTION 



o 

, The possibility to recover its original shape after deformation and heating ["shape- memory effect" (SME)] as well 
as the effect of superelasticity (SE), where elastic strains of up to 8% can be achieved^, make shape memory alloys 
(SMA) to widely used functional materials in industry, ranging from dental to aerospace applications. The physical 
origin of SME and SE is a diffusionless structural phase transition between a high-temperature phase ( "austenite" ) 
| ; , and a low-temperature phase ( "martensite" ) with lower symmetry. In order to explore this kind of phase transitions 
C3 ' theoretically at the nanoscale, molecular dynamics simulations have been performed in the past for different alloys 
£h \ such as NiAl by Saitoh and Lh*£ or FeNi by Entcl and coworkers^. 

I . Among all, the most commonly used SMA's in industry are nickel-titanium alloys with equi- or near-equiatomic 
composition. In these alloys, the experimentally observed ground state structure consists of monoclinic lattice cells 
with additional in-plane and out-of-plane shuffles of the atoms at the faces of the ceh£. This ordered, bi-atomic alloy 
structure is referred to as B19' and has space group P2\/m. If the crystal is heated above a strongly concentration 
dependent tcmperaturo^, it undergoes a structural change to the more symmetric, body-centered cubic B2-phase. In 
addition, the so-called R-structure, an intermediate phase with trigonal symmetry, appears during the martensitic 
B2— >B19'-transition in some experiments^, but there is no overall agreement among them concerning the correct 
space group. By means of ab initio calculations, the phase energetics and lattice parameters of the proposed structures 
have been analyzed in many publications^—. Among those, Huang et al^- were the first, who proposed the body- 
. centered orthorhombic B33 to be the most stable martensitic ground state of NiTi, but this was not yet observed in 
■"sj" ' experiments. 

This ambiguity of the exact crystal structures on the one hand and discrepancies in the predictions of elastic 
constants on the other handi£ make it difficult to construct an interatomic potential suitable for dealing with all the 
observed phenomena in NiTi, especially the structural transitions with the associated shape memory and superclastic 
behavior, by means of molecular dynamics (MD) simulations. That is why only a few potentials exist in the literature, 
but of those, each is capable of describing some aspects of the material properly: Farkas et aL— proposed a potential 
based on the embedded-atom method (EAM), which is able to reproduce lattice parameters and cohesive energies 
of the compounds B2 NiTi and NiaTi, but fails in stabilizing a monoclinic structure in the ground state. With a 
potential originating from the second- moment approximation of the tight-binding model (TBSMA), Lai and Liu2^ 
studied crystallinc-to- amorphous transitions of nickel/titanium solid solutions as well as the amorphization of several 
NiTi compounds upon ion irradiation. Since this potential predicts a monoclinic structure to have a higher cohesive 
energy than the cubic state, it was adopted by Sato et al£^ to perform MD on the stress induced martensitic 
phase transformation, whereby multiple B2-B19' pathways could be identified. Until today, there is a lack of studies 
regarding the temperature driven structural changes in NiTi with simulations. Only Ishida and Hiwatari^ report on 
simulations using the modified EAM, where a reversible phase transition at a certain temperature occurs, but the 
crystallographies of the parent and martensitic phases are not satisfactorily clarified. 

In the present work, the applicability of the TBSMA potential for performing reliable MD simulations of this kind 
of transformation in NiTi is analyzed. For that purpose, the potential is extended by a function, which controls the 
cutoff-behavior of the rapidly decreasing exponential functions appearing in this model, which describe the hopping 
integrals and the pair interaction. This is explained in Sec. Ill Al together with some simulation details. In Sec. Ill Bl 
a method for calculating the free energy in a MD simulation is described, by which the thermodynamics of the phase 



2 



transition can be studied. The modified TBSMA potential leads to a better agreement of lattice parameters and 
energetics of the stable B19'-structure with experiments and ab initio calculations, which is covered in Sec. lIII A"l 
In addition, MD simulations of cquiatomic NiTi, performed under periodic boundary conditions and at varying 
temperatures are presented there, in which hrst order structural phase transformations can be identified, and the 
emerging structures are analyzed. In Sec. lIII 51 the dependence of transition temperatures on Ni/Ti concentration is 
considered and compared with experiments. The strong decrease of these temperatures when the Ni content differs 
slightly from 50% is shown to be attended by a considerable destabilization of the lattice structure. Thermodynamic 
free energy and entropy calculations during the heating and cooling process of the system are presented in Sec. lIII Cl 
In addition, free energy barriers and their dependence on temperature and concentration are calculated along a linear 
transformation path. Finally, a summary and conclusion are given in Sec. lIVI 

II. THEORY 

A. Interatomic potential and simulation method 

The semi-empirical potential applied in this work to perform MD on the martensitic phase transitions in NiTi 
originates from the tight-binding-bond model in second-moment approximation as described by Cleri and Rosato^. 
According to this, the total energy U t of a system is written as a sum over all atoms i: 

u t = Y.(u? + u«), (i) 

i 

with 



and 

The quantum mechanical many-body character of metallic bonding is expressed in the bond energy Iff, which is 
approximated as the negative square root of the second moment of the electron density of states at atom i. This is 
given as a sum over squared two-center- ("hopping"-) integrals, with an effective part and an exponential distance 
dependence. [//* represents a pairwisc Born-Mayer repulsion^, which is necessary to stabilize the crystal, nj denotes 
the distance between atoms i and j, and the indices a/3 account for different atomic types, what leads to 15 parameters 
in binary alloys AB (£, q, A, p, d for A- A, B-B and A-B interaction). 

In the case of NiTi, they were determined by Lai and Liu^S by calculating properties like cohesive energies, lattice 
parameters and elastic constants analytically from the potential and fitting them to ab initio data of the pure materials 
and of the B2-NiTi phase at T = K. Additionally, they treated the cutoff radius of the exponential functions as an 
adjustable quantity, and proposed r c = 4.2 A as an optimal choice. With this value, first and second nearest neighbors 
in NiTi are taken into account in the calculation, but since it lies closely below the distance of the third nearest 
neighbor shell, a smooth cutoff behavior has to be ensured in order to perform reliable MD simulations without 
diverging forces. Therefore, the cutoff function f c proposed by Baskes et al£±: 

{1, r < r c 
0, r>r c 

is multiplied to the exponential functions in the present work, with x = (r c — 5/2 — r) / (5/2). Good agreement with 
structural properties of the monoclinic B19' ground state of NiTi is achieved by setting the additional parameter 5 
to 0.2 A. 

With this potential, MD simulations are performed with a velocity- Verlet algorithm^ (timestep Ai = 10~ 15 s) at 
different temperatures, imposed by a Nose-Hoover thermostat^, and under periodic boundary conditions. A fully 
flexible simulation box is applied by the Parrincllo- Raman method^ with improvements of Martyna et al£^, and the 
pressure is set to Pa. 
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B. Free energy calculation 



The calculation of thermodynamic quantities in MD simulations can give deeper insight into phase behavior and 
phase transitions of many particle systems. While the internal energy is calculated straightforward by taking averages 
of the total energy U t (Eq.[T]), the free energy F has to be obtained in a more complex way, e.g. by thermodynamic 
integration between a reference system with known F and the system under consideration. In the present work, the 
method described by Frenkel and Ladd^S is applied, which makes use of an Einstein crystal (EC) as reference state, i.e. 
a system, where particles are bound to fixed lattice sites by harmonic springs. In the past, this method has successfully 
been used in combination with model potentials to describe e.g. bcc-fcc transitions in iron^i, liquid-crystal interfaces 
in Al22 or phase diagrams of Au-Ni alloys^. 

In detail, a parameter dependent potential U (A) is introduced, with U (1) = U, the potential of the system with 
unknown free energy, and U (0) = Ue, the potential of the EC: 



It can easily be shown, that 



U (A) = U E + A (U - U E ) ■ (2) 
dF(X) ,0Z7(A), 



A 



with the canonical ensemble average (. . . )\. With Eq.[21 this leads to 

F(A=1) = F(A = 0)+ / dX{U-U E )x- (3) 
Jo 

Since F (A = 0), the free energy of the EC, is analytically known for a given temperature T, F (A = 1) can be calculated 
in the simulations. Optimal spring constants for the EC are obtained by averaging and comparing the mean-squared 
displacements for A = and A= 1. In the simulations presented here, this is done at first for 50000 timesteps at a 
given T, followed by a variation of A according to 

, 1 -0.05 _ 1 -0.005 _ M -0.0005 „ 

A : 1 > 0.1 > 0.01 > 0, 

since the integrand of Eq. [3] increases strongly when A approaches 0. At each A-stcp, (U — Ue) is averaged over 10000 
timesteps, and at the end, a numerical integration yields the value of F. As described above, simulations are carried 
out at zero pressure, and therefore F is equal to the Gibbs free energy G, which is used in the following, instead. 



III. RESULTS AND DISCUSSION 



A. Structures and Transformations 



In order to perform simulations of the temperature-driven structural phase transition between B19' [Fig.[TJa)] and 
B2 [Fig.[IJb)] in NiTi, it has to be ensured at first, that the martensitic structure is stable at low temperatures. For 
this purpose, a system with 2048 particles and equiatomic composition was set up in a monoclinic structure with 
lattice parameters a, 6, c and a close to experimental values. Decreasing the temperature stepwise from T = 10 K to 
T = K with a rate of 1 K per 10000 timesteps lead to a stable B19'-structure with parameters listed in Table 1. 

The comparison with experiments and ab initio results shows deviations of about 2%-5% in the a, b and c 
parameters, and an excellent agreement concerning the value of the monoclinic angle a, which is due to the 
modification of the TBSMA potential described above. The in-plane and out-of-plane shuffles with values between 
2% and 9% of the cell parameters b and c do not appear within this model. A calculation of the cohesive energy 
leads to E c (B19') =5.076 eV/atom. 

To check, whether the stable B19' phase undergoes a structural phase transition to the B2 phase upon heating 
above a certain temperature, the system was heated about 1 K every 2000 simulation steps from 2 K to 400 K. In 
Fig. [21 the behavior of the simulation box is recorded during this process. 

The initially set angle a between the a- and c-edge of the simulation box [with respect to the cell shown in 
Fig.QJa)] first increases due to thermal expansion, but then drops immediately to 90° when reaching 318 K, whereas 
the angle /3 between the a- and c-edge stays nearly constant. At the same time, the different lengths of the a- and 
6-cdges of the box reach the same value, which lies almost at the arithmetic mean. Since in addition neither the 
angle between the b- and c-edge (90°) nor the length of the c-edge changes, a cubic box is obtained after the phase 
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FIG. 1: (a) The monoclinic B19' structure of NiTi with angle a, lattice constants a, b, c, and shuffles, denoted by the arrows; 
(b) the cubic B2 structure. 



TABLE I: Calculated values of structural parameters (see Fig.Q]) and the cohesive energy difference to the B2 structure in 
comparison with results for the unmodified TBSMA potential, ab initio calculations and experimental values. The lattice 
constants a, b, c are given in A, the energies in meV/atom. 



Structure a b c a AE C 

B2 (this work) 3.01 90.0° 

B2 a 3.01 90.0° 

B2 b 3.019 90.0° 

B2 C 3.013 90.0° 

B19' (this work) 4.45 4.03 3.00 97.7° 54.0 

B19' a 4.46 4.19 2.96 93.3° 18.0 

B19' fc 4.677 4.077 2.917 98.0° 55.42 

B19' d 4.66 4.11 2.91 98.0° 



"Calculated results with the unmodified tight-binding potential by Lai and Liu-*. 
^Calculated ab initio results by Hatcher et a/.—. 
Experimental results by Sittner et aZ.— . 
Experimental results by Prokoshkin et oZ.— . 
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FIG. 2: Shape changes of the simulation box during a temperature cycle: (above) angles a (between a- and c-edge) and j3 
(between a- and 6-edge) as well as the imposed temperature. (Below) lengths of a- and 6-edges . 



transition. Assuming a B2 structure, the lattice parameter c can be calculated from the box dimensions, yielding 
c = 3.01 A. With this value, the cohesive energy of a perfect B2 lattice is computed within the model potential 
(Eq.HJ to E c (B2) = 5.022 eV/atom, which leads to E c (B19') - E c (B2) = 54 meV/atom in good agreement with ab 
initio results^ (see Tab.|l|. As the system is cooled down again, the box transforms with a hysteresis of 28 K to a 
shape, which differs from the starting B19' geometry. 

To get more insight into the crystallography of the involved structures, the radial distribution function g (r) is 
evaluated at different points in the temperature cycle, as shown in Fig. [3] In the heating process at 100 K, where the 
phase transition has not occurred yet, the peaks of g (r) lie at the positions of the perfect B19' starting configuration, 
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with the typical broadening due to thermal vibrations. Above the transition temperature, at 400 K, the structural 
change can be identified by the appearance of a peak at r = 4.3 A and a lowering of the peaks at r = 4.47 A and 
r = 4.75 A. The comparison of this structure with the previously assumed perfect B2 shows good agreement, only 
the peak at r = 4.04 A indicates a slight difference. The transition occurring in the cooling process finally leads to a 
structure, for which the radial distribution function g (r) shows nearly the same peaks as in the starting B19' lattice, 
but since g (r) cannot resolve the angular distribution of the appearing lengths, the identity of the structures cannot 
be concluded. 
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FIG. 3: Radial distribution function at different temperatures in the heating and cooling process together with the positions 
of the perfect structures. 



In order to clarify the crystallographic details, the nearest neighbor environment of the low- and high-temperature 
structures is analyzed in detail by averaging the neighbor distances of each atom over 10000 simulation steps. This 
is done at low temperatures before heating (B19') and after cooling (B19'), and at T = 350 K (B2). The results for 
two characteristic atoms with nearest neighbors are shown in Fig. [4] 

In the stable starting configuration B19', two of the eight nearest neighbor (NN) distances are elongated about 
13.5% [dashed lines in Fig|4ja)] relative to the others due to the monoclinic shear along [100]. In an ideal B2 structure 
with c = 3.01 A all the NN lengths would have the same value of 2.61 A. In fact, the high-temperature phase obtained 
in the simulation, shown in Fig.@|b), consists of 4 NN lengths of about 2.54 A and 4 NN lengths of about 2.69 A. 
These deviations of ± 3% of the perfect B2 value cause the slight differences in the radial distribution function 
(Fig. [3]), but since they occur in an alternating manner, the simulation box is nevertheless cubic and the structure 
is closely related to a perfect B2. If the system is cooled down again, the emerging structure below the transition 
temperature [Fig.HJc)] consists of the same, but resorted NN distances as in B19' [Fig.HJa)], which explains the 
nearly identical peak structure of the radial distribution function but a change in the shape of the simulation box. 
This structure has a cohesive energy of E c = 5.078 eV/atom, which is only about 0.04% higher than E c (B19'), and 

it will be denoted as B19' in the following. Therefore it can be stated, that B19' is very close to the global cohesive 
energy maximum within the used model. It is the authors opinion, that the small deviations of structures and lattice 
constants from experimental results and ah initio calculations have to be traced back to the TBSMA potential model, 
which is an approximation and simplification of the real force field in NiTi. 

In contrast to experiments, where a strong twinning of different martensitic variants is observed in NiTi nanocrys- 
tals when cooled down below the martensitic transition temperature^, the ground state of the simulated system 
(B19') consists of one single variant. This may be due to the periodic boundary conditions and the simulation box 
as well as a system size, which is much smaller (« 3 — 4 A) compared to the nanocrystals (w 50 A). 

Since all the simulations are performed in a box with 2048 particles and periodic boundary conditions, it 
has to be ensured, that finite size effects do not play a decisive role concerning the structures, energetics and 
transition temperatures. To this end, the system size was varied between 500 and 32000 particles, whereupon 
the only differences were observed in shifts of the transition temperatures Ta (B19'— >B2) of about 5 K and 
Tm (B2— >B19') of about 10 K when going to the larger systems (see Fig. [5]). Furthermore, the involved struc- 
tures and structural energies do not show a dependence on system size. The hysteresis between austenitic and 
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FIG. 4: Nearest neighbor (NN) environment of the emerging structures: (a) B19' with NN distances of 2.51 A (continuous) 
and 2.85 A (dashed); (b) high-temperature B2 like phase with NN distances of 2.54 A (continuous) and 2.69 A (dashed); (c) 
ground state structure emerging upon cooling of the B2 like lattice, lengths as in (a). 

martcnsitic transition takes a constant value when the system size exceeds 7000 particles, Ta — ?m = 35 ± 5 K, 
in good agreement with recent experimental values^ (see also Scc. lIIIBj) . Therefore the existence of free energy 
barriers and associated hysteresis between low- and high-temperature structures is not an effect of finite system 
size, but rather a principal attribute of first-order martensitic phase transitions (see Sec. lIII C|) . So it is reason- 
able to compare the transition temperatures resulting from the simulations with experimental values in the following. 
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FIG. 5: Dependence of austenitic (Ta) and martensitic (Tm) transition temperatures on system size. The error bars result 
from a temperature step of 5 K during heating and cooling. 



B. Dependence on concentration 



It is a well-known experimental fact, that temperatures, at which the structural changes between martensite and 
austenite in NiTi take place, are very sensitive to small deviations of the Ni/Ti concentration from the ideal one with 
50% Ni and 50% Ti£i2Zr— . The effect, that B19'— S-B2 transition temperatures in NiTi decrease with increasing Ni- or 
Ti-conccntration is reasonable in a sense, since the crystallographic ground states of the pure materials Ni (FCC) and 
Ti (HCP) differ from a monoclinic structure corresponding to the alloy phase B19'. If the concentration of one of the 
constituents is rised, more and greater islands of pure material emerge, which are fixed in a matrix of an energetically 
less favored structure. 

In order to confirm this behavior theoretically within the semi-empirical approach, simulations were carried out 
where a temperature cycle was applied to systems with 2048 particles (as in Sec. lIII A]) and different Ni concentrations 
between 47% and 53%. Starting from a perfect NisoTiso in the B19' structure, these compositions were achieved 
by substituting a commensurate amount of atoms of one sort by atoms of the other sort randomly. In the heating 
process, the B19'— ^B2 transition occurs instantaneously at a temperature denoted by Ta, and upon cooling, the 
system transforms into the B19' structure at the temperature Tm- For each concentration, 10 differently assembled 
systems were simulated. The averaged results are shown in Fig.|6l together with experimental results for comparison. 
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In the experiments, a difference is observed between the temperatures where the martensite nucleation starts (M s ) 
and ends {Mf), and where the austenitic transition begins (A s ) and ends (Af), which leads to the definitions of 
T M = 0.5 {M s + M f ) and Ta = 0.5 (A s + A f ). 
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FIG. 6: Dependence of the transition temperatures on Ni concentration. Simulated data (a) in comparison with experimental 
values (bjaa, (c)4, (d)2&, (e)^. 

On the Ni-rich side, transition temperatures agree well with recent results of Khalil-Allafi and Amin-Ahmadi^ 9 -. 
Since they measured differences M s — Mf and Af — A a in the range of 29 K to 45 K, the simulated curves lie nearly 
in between these values. In addition, the extent of hysteresis, Ta — Tm ~ 30 K, coincides well. For Ni concentrations 
c N i > 51%, the measured decrease of M s seems to be yet steeper than in the simulations, which is confirmed by 
the experiments of Hanlon et al&. It was pointed out by several author o 35 ^ 38 ' 40 ! 41 , that experimental transition 
temperatures in NiTi depend strongly on the processes of material preparation before the actual measurement, since 
varying ageing temperatures and durations as well as quenching rates can lead to precipitation effects of intermediate 
phases. Experiments of Saburi et at— for an alloy with c Ni = 51.3% show for example M s values between 240 K and 
77 K, depending on heat-treatment. 

The discrepancies between the simulated transition temperatures and the measured values at the Ti-rich side, where 
only very few experiments exist in the literature, are possibly due to Ti-precipitation effects in the measured samples, 
since lattice sites, where Ni atoms have been replaced by titanium, are attractive to each other—. Therefore, they 
can form larger Ti islands influencing the transition behavior and temperatures. In the simulations, the replaced Ti 
or Ni atoms stay at their positions, leading to nearly symmetric TT curves around 50% nickel. 

In order to explore the reasons for the strong dependence of Ta and Tm on c Ni , the local environment of the atoms 
is examined. As shown in Fig-[4fa) and (b) (see Sec. MI A"| . a structural feature for distinguishing unambiguously 
between B19' and B2 locally is the mean length of the two nearest neighbor distances in [101] and [101] direction, 
since the elongations due to the monoclinic shear in B19' vanish when the system transforms to B2. Defining this 
length as d leads to d (B19') = 2.85 A at T = K and d (B2) = 2.61 A at T > T A in Ni 50 Ti 5 o within the used potential 
model. Information about the local structures of all atoms in the system can be gained by evaluation of n(r), the 
number of atoms with d~r divided by the number of all atoms. This quantity is shown in Fig. [7] at four different 
temperatures during a heating process with starting configuration B19', for systems with Ni concentrations of 50%, 
51% and 52%. 

At T = 10 K, the system with 50% Ni shows a sharp peak at r = 2.85 A, as expected, whereas one additional percent 
of Ni lowers this height by a factor of about 0.5, what is proceeded when going to 52% of Ni. Because of a constant 
number of particles, a lowering in height is equivalent to a broadening of the peaks, or more and more lengths d 
differing from the ideal one at T = K. This is located at the peak maximum, which itself is shifted by about 0.2 A 
per Ni-percent, resulting in smaller monoclinic shear angles a. Heating the system leads to thermal peak broadening 
and small horizontal shift, but the effect of concentration to the relative peak structure remains unchanged. So it can 
be stated, that a small deviation from c Ni = 50% is attended by a strong destabilization of the perfect B19' lattice 
structure, independent of temperature, which can be regarded as a reason for the strong decrease of Ta with c Ni . 
At T = Ta, all atoms of the system transform collectively to B2 without any preceding structural changes, since the 
peaks are still located around d«2.85 A at temperatures slightly below Ta, but shifted completely to <i«2.61 A at 
T > Ta, after the phase transition. 
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FIG. 7: Local structure analysis via the number of atoms with d = r divided by the number of all atoms, n (r), for 50%, 51% 
and 52% Ni, at 3 temperatures below Ta and one greater than Ta- 



Further insight into this destabilization effect can be gained by calculating the equilibrium cohesive energies for 
different Ni percentages between 48% and 52%. For this purpose, the system was set up in B19' structure and the 
temperature was reduced stepwise to K. The results are shown in Fig. El together with the contributions of bond 
energy and pairwise interaction (Eq.QJ. The energy references are set to the values of the perfect structure with 
50% Ni. Varying the Ni percentage results in reduction of the cohesive energy in a symmetric way, which is strongly 
related to the decrease of transition temperatures (Fig. [5]). This behavior is due to the bonding part of the energy, 
which explains more instable lattices at Ni concentrations varying from 50%. 
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FIG. 8: Dependence of the equilibrium cohesive energy E c , bonding part and pairwise interaction part of E c on Ni concentration. 
Energies are shown relative to the system with 50 % Ni. 

Considering a perfect structure (p) on the one hand, and a structure with an "impurity" (i), where one Ti atom is 
replaced by a Ni atom on the other hand, cohesive energy differences AE C = EP — E l c of individual atoms at and around 
the impurity can be calculated. For the impurity atom itself, AE C = 0.523 eV is obtained, as well as AE C — 0.026 cV, 
if the impurity is a nearest neighbor atom, and AE C = 0.006 cV in the case of a second nearest neighbor. This shows, 
that there is a lower value of E c in the direct vicinity of an additional Ni atom as in the surrounding perfect B19' 
NiTi, which destabilizes the structure locally. 

Further inspection of the bond energy, according to the tight-binding parametrization (see Eq.[T]), reveals, that by 
replacing a Ti atom by a Ni atom, the parameters q, £, and d change in such a way, that the bonding part of the 
cohesive energy of this Ni- "impurity" is lower than for the Ti atom. Thus, the "impurity" scenario is energetically less 
favorable, resulting in an effective "repulsion" compared to the impurity-free case, destabilizing the lattice. Similar 
arguments apply for the Ti-rich side of the phase diagram, which explains the decrease of transition temperatures in 
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this part of the phase diagram. 



C. Thermodynamics 



Martcnsite/austenite phase transitions can in general be explained by means of different competing contributions 
to G, the Gibbs free energy. One phase can only form within a matrix of the other one, if non-chemical clastic strain 
energies and energy differences due to unlike interfaces (AG nc ) are overcome by a chemical energy AG^. Therefore, 
the structural transition occurs, if AG C = AG nc , and leads to a jump of G about this amount, which characterizes 
the transition as first-order. Thermodynamic properties of the phase transitions visible in the presented simulations 
are calculated by thermodynamic integration as described in Sec. Ill Bl A 2048 particle system in B19' configuration 
is heated above Ta to T = 360 K and subsequently cooled. At temperature intervals of 20 K (10 K in the vicinity of 
the structural change), the Gibbs free energy is determined. Together with the internal energy U and the entropy 
S ={U -G) /T, G is shown in Fig.® between T= 150 K and T = 360 K. 
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FIG. 9: Energy U, Gibbs free energy G and entropy S between T = 150 K and T = 360 K during a temperature cycle. 



At T = Ta, G drops about AG C ~9 meV/atom from Gbw to the lower Gb2- Since at this temperature U increases 
stronger than only due to thermal effects as in the previous heating, the negative AG has to be accompanied by a 
jump in entropy to a higher value. Khalil- Allan and Amin-Ahmadi measured austenitic start- and end-temperatures 
(A s , Af) as well as enthalpies AHm->a hi B19'— >B2 transitions for NiTi with Ni content between 50.3% and 51.0%2£. 
Ta = 0.5 (A s + Af) and AHm-^a show nearly linear behavior with composition in this range and can therefore be 
extrapolated to 50% Ni, what leads to ASm->a = AHm-^a/Ta = 0.044 meV/(atom-K), in good agreement with the 
value obtained in the simulation of 0.045 meV/(atom-K). Several ab initio studie a 12 i 42 ~ — state an entropic stabilization 
of the high temperature phases in NiTi by vibrational entropy and phonon mode softening, respectively. If, in the 
present simulation, the temperature is decreased again, a jump of G at T = Tm back to the B19'-curve is not observed, 
because the system does not transform to B19' upon cooling, as explained in Sec. MI Al Nevertheless, investigation 
of U and S shows the occurence of a phase transition (B2— >B19') by revealing the contrary behavior to the heating 
process: the energetically more favorable structure with lower entropy is adopted by the system. 

Moreover, with the help of free energy calculations it is possible to detect free energy barriers in transformation paths 
between the occurring structures, which are typical in first-order phase transitions and responsible for hysteresis. As 
an example, a linear path between B19' and B2 is modeled, and at each intermediate step, the free energy is calculated 
by a simulation as above. To this end, the system is set up in a structure between B19' and B2 with edge lengths 
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a, b, c and an angle a of the unit cells [see Fig.flja)] according to: 

o(e) = a B i9' + e • (a B2 - aBi9') 
b{e) = 6bi9' + e • (&B2 - &B19') 

c(e) = C B 19' + £ • (c B 2 - C B 19') 

ot(e) = a B i9' + e • («B2 - «Bi9') 

with a parameter e varying from (B19') to 1 (B2), and values for the reference structures obtained in this work (see 
Tab. |l| . The angles j3 (between a- and &-edge) and 7 (between b- and c-cdge) remain at 90°. The system is forced to 
stay in this structure by a stiff simulation box with edge lengths and angles according to the unit cells values, but the 
larger the system size, the less is a stiff box able to hold the system into a defined structure. Therefore the number of 
particles is reduced to 500 in this study. Fig. ITUl shows the results at different temperatures T and Ni concentrations 
c Ni in form of AGy. CNi (e) = Gj\ CNi (e) — Gt iCN1 (0)- Thus, the differences in G can be seen along the transformation 
path and compared for varying T and/or c Ni , without accounting for the absolute G values. 
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FIG. 10: AGr,c Ni (e) = Gt, cn1 (e) — Gr,c Ni (0) along a linear transformation path from B19' (e = 0) to B2 (e = 1) for temperatures 
T=100 K, 200 K, 300 K and nickel concentrations c Ni = 50%, 51%. The error bars at the values for c Ni = 51% result from 
averaging over 3 randomly assembled systems. 



Regarding the system with c Ni = 50%, a free energy barrier can be detected for T=100 K, which flattens when 
going to higher temperatures. The B2 structure lies in a minimum of AG, too, but this is not surprisingly, since 
the potential parameters were determined by fitting to material properties of B2 at T = K. Nevertheless, at low 
temperatures, B19' exhibits a lower AG, but already at T = 200 K, more than 100 K below the phase transition, B2 
would be more stable, and only the barrier in AG prevents the system from transforming. Upon cooling, the barrier 
does not disappear, which explains, why a B2— »B19' transition is not observed in the simulations. Increasing the Ni 
content to 51% leads to a less pronounced free energy barrier at T= 100 K with a more kinked curve shape, which is 
caused by slight reorientations of the initially set up structures, and to considerable uncertainties when averaging over 
3 differently assembled systems. This behavior reflects the result of Sec. lIIIBl where a more unstable and thus more 
fluctuating lattice structure has been obtained by varying the Ni concentration away from 50%. As a consequence, 
transition temperatures decrease, what is confirmed here, too, since the AG barrier vanishes already at about 200 K. 



IV. SUMMARY AND CONCLUSION 



In this work, MD simulations of the temperature driven structural phase transitions in NiTi alloys with equi- or 
nearly equiatomic composition were carried out by using a semi-empirical model potential from the literature^ (TB- 
SMA), since this is known to predict a monoclinic structure to be energetically more favorable than a cubic one at 
T = K. 

It could be shown, that a B19' structure is stable within this model at low temperatures, and a slight modification 
concerning the cutoff behavior of involved functions leads to good agreement of lattice parameters and energetics 
with ab initio and experimental results. By analyzing simulation box shape, radial distribution function (RDF) and 
nearest neighbor (NN) environments during an increase followed by a decrease of the imposed temperature, structural 
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phase transitions were observed. While at T = Ta during heating (318 K for N150T150) the system adopts a nearly 
cubic structure closely related to B2, at a temperature T — Tm upon cooling, a structure denoted B19' emerges. This 
has nearly the same RDF than B19', but the NN's belonging to the RDF peaks are resorted, resulting in a larger 
cohesive energy of about 0.04%. So it can be stated, that the TBSMA approach predicts a stable B19', which is very 
close to the martensitic ground state B19' of this model. 

The experimentally known fact, that transition temperatures (TT) vary strongly with nickel concentration was 
confirmed qualitatively at the Ni-rich side and even quantitatively in the range between 50% and 51% nickel. At the 
Ti-rich side, there exist only a few and partly conflicting experimental results, which do not fit well with the simulated 
curves. The discrepancies may be due to precipitation effects in the material processing before the measurements, 
which do not emerge in the simulations. By investigation of NN distances for systems with c Ni > 50% during heating, 
it could be shown, that the strong decrease of the austenitic TT is attended by a dcstabilization of the B19' lattice 
structure. This destabilization results from the bonding part of the cohesive energies in the direct vicinity of a Ni- 
"impurity", which is lower than in a perfect B19' NiTi. A relation between lattice stability and TT values was also 
proposed by Lu et al*^, who performed ab initio charge density calculations of B2 Ni-rich NiTi. 

Thermodynamic calculations gave more insight into the phase transition process by confirming that the high tem- 
perature phase is entropically stabilized, with a jump of AS at T = Ta fitting well with recent experimental results^. 
Free energy barriers suppressing the phase transitions were determined along a linear path between B19' and B2, and 
it could be shown, that the heights of these barriers are lowered by increasing on the one hand temperature and on 
the other hand nickel concentration from 50% to 51%. 
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